Two-dimensional Particle-in-Cell modeling of blow-off impulse by X-ray irradiation* 
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Space objects, such as spacecraft or missiles, might be exposed to intense X-rays in outer space, leading to 
severe damage. How to reinforce these objects to reduce damage from X-ray irradiation is a significant concern. 
Blow-off impulse (BOT) is a crucial physical quantity for investigating the material damage induced by X-ray 
irradiation. However, the accurate calculation of the BOI is a challenge, particularly for the large deformation 
of materials with complex configurations. In this paper, we develop a novel two-dimensional Particle-in-Cell 
(PIC) code, Xablation2D, to calculate the BOIs under far-field X-ray irradiation. This significantly reduces the 
dependence on grid shape for numerical simulation. The reliability of this code is verified by the simulation 
results from the open-source codes, and the calculated BOIs are consistent with experimental and analytical 


results. 
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I. INTRODUCTION 


High energy-density X-ray irradiation induced by a nuclear 
explosion in outer space can damage spacecrafts or missiles 
[1-3]. If the X-ray source is located in the far-field, the sur- 
face material of these objects might undergo sublimation. The 
vaporized material rapidly expands outward, causing a blow- 
off impulse (BOI). The BOT loads on the remaining solid ma- 
terial, generating a compressive stress wave that propagates 
inward, which is called vapor recoil loading. On the other 
hand, the pressure disturbance caused by the non-uniform de- 
position energy as thermal stress loading generates a thermal 
shock wave characterized by compression at the front and ten- 
sion at the rear [4]. These stress waves form a thermal shock 
wave, and induce mechanical damage to material structures 
and make space objects lose efficacy permanently. How to 
reinforce the space objects to reduce damage by X-ray irradi- 
ation is a valuable issue that has been investigated extensively. 

The BOI, as a measurable physical quantity, plays a crucial 
role in investigating the material damage induced by X-ray 
irradiation [5—8]. Since the pulse duration of X-rays induced 
by a nuclear explosion and the time taken for phase transition 
in the material are both O(0.1 us), significantly shorter than 
the time it takes for dynamic response and stress wave propa- 
gation within the material. It is reasonable to decouple the en- 
ergy deposition process from the whole X-ray irradiation. In 
early research, several physical analytical models were pro- 
posed to calculate the BOI under the assumption of instan- 
taneous deposition energy [9, 10], including the Whitener 
model, the BBAY model, and the modified BBAY (MBBAY) 
model. These models have analytical formulas for the BOI, 
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and they are extensively utilized in subsequent simulation and 
experimental works [9, 11-14]. The predictive capabilities 
of these models are satisfactory when it comes to calculating 
the BOI for one-dimensional or simple configuration mate- 
rials, but their accuracy falls short when applied to complex 
configuration materials and non-instantaneous deposition en- 
ergy. The other approach for calculating the BOI is to com- 
bine the energy deposition process with the generation and 
propagation of stress wave, and simulate the whole X-ray ir- 
radiation. With the rapid development of computational fluid 
dynamics, a series of codes, referred to this approach, have 
been developed to calculate the BOI for multi-dimensional 
materials, such as PUFF-TFT [15], CTH [16], LS-DYNA [17], 


ss ABAQUS [18], RAMA [19-21], TSHOCK3D [22] and so on. 
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These codes employ either the finite difference or finite el- 
ement method, and incorporate Eulerian or Lagrangian de- 
scriptions for numerical simulation, which enables the cal- 
culation of the temporal evolution of the BOI and results in 
significant enhancements in computational accuracy. Never- 
theless, there are certain disadvantages associated with these 
codes, i.e., the Eulerian description is hard to track fluid in- 
terfaces, while the Lagrangian description is prone to mesh 
distortion when encountering large material deformations. In 
addition to the mesh method mentioned above, another simu- 
lation tool is the Monte Carlo method [23, 24]. However, this 
method relies on statistics, has low efficiency, and requires a 
significant amount of computational resources. 


Particle-in-Cell (PIC) is a particle-mesh method that is de- 
veloped by F. H. Harlow’s team for the first time when study- 
ing gas dynamics problems at the Los Alamos National Labo- 
ratory in the United States [25, 26]. Then, it is widely gener- 
alized to computational plasma physics [27]. The Harlow’s 
PIC method discretizes the fluid into free-moving pseudo- 
particles in the spatial grid, and combines the Lagrangian 
and Eulerian descriptions, which has advantages in simulat- 
ing large deformation problems in materials with complex 
configurations. As far as we know, there is no reports on 
PIC code for the BOI calculation. In this paper, we develop 
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a novel two-dimensional PIC code, Xablation2D. This code 
can be used to calculate the BOI produced by X-ray irradiat- 
ing materials. 

The paper is organized as follows. In section II, we dis- 
cuss the theoretical basis and essential modeling techniques, 
including the fluid PIC scheme, and the three main modules 
of the code. In section II, we describe the algorithm im- 
plementation in detail, including initialization, discretization, 
and parallelization. Section IV is dedicated to present the 
reliability of the code, which includes the simulation of the 
energy deposition, the shear flow and the Kelvin-Helmholtz 
instability. Furthermore, a series of material simulations in- 
volving X-ray irradiation are conducted in this section to ver- 
ify the simulation capabilities of the Xablation2D code for the 
BOI. Finally, we arrive at a conclusion in the last section. 


Il. THEORETICAL BASIS AND MODELING FOR 


SIMULATION 
A. Fluid PIC scheme 


We use the Harlow’s splitting PIC scheme to solve the hy- 
drodynamic equations. Without loss of generality, the equa- 
tions can be written as an abstract operator form 


() 


where q(r, t) is an arbitrary hydrodynamic quantity, and Ais 
an abstract operator. According to the rule of operator decom- 
position, the solution of the equation (1) at the time step At is 
reduced to the sequential solution of two auxiliary problems 
[28] 


Oa) + Palri) = 0, 
ae t) 2) 
2d + fa(r,t) =0, 


where A = É + Ê. The two equations in Eq. 2 correspond to 
the "Euler step" and "Lagrange step" in the Harlow scheme, 
respectively. In the “Euler step”, the operator Ê does not in- 
corporate the spatial divergence operator, so the equation is 
easy solved on a fixed spatial grid. In the “Lagrange step”, 
pseudo-particles are introduced to carry mass density, mo- 
mentum density, and specific internal energy density. In one 
time step, as the pseudo-particles move to new positions, the 
new hydrodynamic quantities are obtained on grids by sum- 
ming up pseudo-particles. The “Lagrange step” can be seen 
as a computational procedure for modeling particle migration, 
which compensates for the transport effect that is neglected in 
the "Euler step". 

Without loss of generality, the second equation in Eq. 2 can 
be written as 


ôq 


At + V-(qU) =0, 


(3) 
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where U is the flow velocity, and q = (p, pU, pe) repre- 
sents the mass density, momentum density, and internal en- 
ergy density, respectively. Eq. 3 is the form of the conserva- 
tion equation, the solution of which can be written as a sum- 
mation of pseudo-particles 


N 
q(r,t) = D Q;R(r,r;(t)). (4) 


The total number of pseudo-particles is denoted as N, Q; 
represents the carried hydrodynamic quantities of the j-th 
pseudo-particle. The value of Q; remains a constant in one 
“Lagrange step”. The kernel function of the pseudo-particle, 
denoted as R(r,r;(t)), is a function that depends on the cur- 
rent space coordinate r and the radius-vector r; of the j-th 
pseudo-particle center. This kernel function satisfies certain 
universal properties in the usual 


aR __ƏR 
Ory Ora’ (5) 
dr, R(r1, r2) = 1, 
Q 


where Q denotes the full space. With the aforementioned 
characteristics and any smooth finite function, the represen- 
tation q(r,t) in Eq. 4 enables the simplification of Eq. 3 to 
fulfill the equations of motion for pseudo-particles 


dr j 

dt 

Our code comprises three primary modules: energy de- 
position module, equation of state (EOS) module and ideal 
hydrodynamics module, which are complemented by a post- 


processing script for blow-off impulse calculation to consti- 
tute the complete Xablation2D code. 


=U(r;(t)), j=1,2,...,N. (6) 


B. Energy deposition 


The energy deposition module is responsible for estimating 
the transfer of energy between X-rays and matter. At a mi- 
croscopic level, X-rays primarily interact with matter through 
electrons. Initially, the energy of the photons is transferred di- 
rectly to the electrons, which then interact with atoms in the 
matter, resulting in energy deposition. Photoelectric effect 
dominates in the low energy region. With the energy of the 
photons increases, Compton scattering would contribute [29]. 
In the case of far-field X-ray irradiation the degree of ioniza- 
tion of the material is so low that effects of plasma could be 
disregarded. The energy flux for parallel X-rays incident on 
the target material can be estimated by considering the dis- 
tance x traveled in the incident direction, 


Sp Hi f(A, T) exp|-u(à)pz]dà 
I Fen 


P(x) (7) 


where ®o and (x) are energy flux at initial and x position 
respectively, f(A, T) is the energy spectrum for X-ray, which 


155 depends on the photon wavelength À and radiation temper- 
iss ature T, p is the mass density, (A) is the mass absorption 
157 coefficient associated with the photon wavelength. The en- 
ergy spectrum f(A, T) approximates a black-body spectrum 
for X-rays produced by a nuclear explosion. In then numerical 
simulation, it is necessary to truncate and discretize energy 
spectrum. Hence, the expression of energy flux in computing 
can be written as, 
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163 


B(x) = Bo >, wj exp|-u(à;)ox], 


where subscript j is the discretized energy group index, 
wj represents the proportion of incident energy flux of 
monochromatic light with a specific wavelength Aj. The 
quantity er signifies the amount of photon energy deposited 
through the interaction between X-rays and matter, per unit 
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mass and per unit time, within the region x ~ x + Ax 


P(x) — O(a + Az) 
pAxt 


a= (9) 


where 7 is the time increment. It is postulated that the en- 


tirety of the aforementioned deposited energy undergoes con- 
version into internal energy in the subsequent discuss. 


174 C. Equation of state 


X-ray irradiation induces significant changes in the state 
of matter, including phase transitions, thermal expansion, and 
shock compression etc, rendering a large range of material 
178 parameters. Consequently, it is necessary to employ distinct 
179 equations of state to characterize the expansion and compres- 
180 Sion regions of the material, respectively. For the thermal ex- 
181 pansion region, the PUFF EOS [15] is adopted, 


182 


183 Where po and p are the initial and current density respectively, 
184 p is the material pressure, I'o is the Giineisen coefficient, ~y 
is the specific heat ratio of vaporized gas, e, is the sublima- 
tion energy, and N = Cĝ/Toes. Co is a Hugoniot parameter, 
which determines the shock wave velocity D in solid material 
with the post-shock velocity u and another Hugoniot param- 
eter by 
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190 D = Co + àu. (11) 


For the compression zone, the Güneisen EOS on the Hugoniot 
line is used, which is expressed as [30] 
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192 


193 p = pH (v) + polo(e— en), PÈ po, (12) 
14 where py (v) and ep are 
C2(1— v/v 
486 pi(v) = Po 0 ( / o) (13) 
[1 — A(L — v/v0)"] 
196 and 
1 
197 eH = gH (vo = v), (14) 


respectively, which represents the post-shock pressure and 
specific internal energy when the pre-shock is stationary, re- 
spectively, v = 1/pis the specific volume, and vp is the initial 
specific volume. 

In addition to the analytic expression of EOS, an open- 
source code, FEOS, providing the EOS for a wide range of 
temperatures and densities in tabular form [31]. FEOS based 
on the QEOS (quotidian equation of state) model [32], cal- 
culates the material specific Helmholtz free energy F'(p, T) 
207 directly, which is widely used in some computational fluid 
28 dynamics codes. 
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D. Ideal hydrodynamics 


209 


In this hydrodynamics module, on the one hand we 
1 ignore the effects of thermal radiation and heat conduc- 
tion. The far-field X-ray flux we concerned on the issues 
is O(100 J/cm?). In this case the material temperature 
is O(1 eV), and corresponding radiation pressure is only 
O(1 Pa), which is far lower than the material pressure. 
The velocity of the thermal shock wave in solid materials is 
about O(1 km/s), which is much faster than heat conduc- 
tion. Therefore, it is reasonable to ignore these effects. In the 
other hand, our code aims at calculating the blow-off impulse 
occurring in the thermal expansion vaporized region. This re- 
gion is characterized by significantly lower material stress 
compared to the material pressure, thus we can also ignore 
the stress for impulse calculation in the following. 

The 2D governing equations of ideal hydrodynamics are 
written in the following, 
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227 where p is the mass density, u = (uz, wy) is the fluid velocity, 
228 p is the pressure, e is the specific internal energy, w = p(e + 
229 U? /2) is the total energy, and ep is the deposited energy from 
230 X-ray irradiation per unit mass per unit time. Account to the 


231 Harlow’s splitting scheme, the 2D governing equations can 260 
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where subscript 1 and 2 denote “Euler step” and “Lagrange 
step”, respectively. 


E. Blow-off impulse calculation 


When X-ray irradiates on the material surface, the material 
undergoes sublimation, forming an evaporation layer. The va- 
porized material violently ejects outward to the surrounding 
and generates a blow-off impulse. According to the impulse 
theorem, the specific impulse in a given direction equals to 
the change in momentum 


P 
ya 
Po 


where P, u, Po and uo represent the final momentum, final 
flow velocity, initial momentum and initial flow velocity in a 
certain direction respectively, m is the material mass. Typi- 
cally, the initial velocity of the vaporized material is zero, that 
is uo = 0. For numerical calculations, Eq. 18 can be written 
as a summation form [33] 

2 


uj<0,ej>es 


dP = P — P) = m(u — uo), (18) 


I= (19) 


m,|ujl, 


where my and u, are the material mass and velocity in the j-th 
grid respectively, that accounts for the sum of the momentum 
in all grids in which vaporized material ejects outward. On 
the other hand, the BOI can also be computed by its definition 


t 
i=] pgdt, 
to 


where p, represents the pressure exerted by the ejected gas on 
the surface of the solid material at rest, t and to are the final 
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and initial time, respectively. For numerical calculations, the 
discrete form of Eq. 20 is 


I= X pAt, 


where n represents the n-th vaporized grid. 


(21) 


II. ALGORITHM IMPLEMENTATION 


According to the theoretical framework outlined in Sec- 
tion II, we develop a two-dimensional code, named as Xab- 
lation2D, to calculate the blow-off impulse of material under 
far-field X-ray radiation. This section provides a detailed de- 
scription of the algorithm implemented in Xablation2D. In 
the following discussion, we will omit the subscript 1,2 that 
distinguishes “Euler step” and “Lagrange step”, for conve- 
nience. 


A. Initialization 


We discrete the simulation domain into rectangular grids 
on the Cartesian coordinate. The physical quantities of the 
fluid, exception of mass density, are initialized by manually 
assigned to the grid points. The mass density is initialized 
by summing up the weight of the pseudo-particles shown as 
follows, 


Pi,g = Pij + Pp? Wig, 


Pi+1,j = Pi+1,j T Pp’? i,j; (22) 


Pi,jtl = Pi,j+1 T Pp Wijs 


Pi+1,j+1 = Pi+1,j+1 F Pp * i,j; 


where p; j is the mass density on the grid, pp denotes the mass 
density carried by the pseudo-particle, and w; j is the weight 
factor. It is important to note that in our code, the value of pp 
remains constant throughout the entire simulation. The area- 
weighting method shown in Fig. 1 is employed to distribute 
the mass density in the 2D simulation, and weight factors are 
expressed as 


“ii = A; + Az + A3 + Ay? 
witii = a= A + A; 1 A’ a 

“STE A; F Az + Aga Aa’ 

MELIH Ar pAg kA FA 


where A; (i = 1,2,3,4) represents the area of overlap be- 
tween the pseudo-particle and the grid. 
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Fig. 1. Area-weighting method for the pseudo-particle in 2D simula- 
tion.The circle p represents the pseudo-particle center, and the black 
dot (i, 7) is the grid point. 


B. Euler step 


The radiative deposition energy and equation of state are 
discretized on the grid in the following, 


n n 
_ Biaj- Pij 


CR ig = ) (24) 


Pi jħaT 
and 
Pij = PHij 
+Topo | e745 — ZP H,i,j (Vo Vig) | > Pij 2 Po 
Pij = Pij 


where ph ig and e"4 i,j are 


no poCg(1 — vi; /vo) 
[1 — A(1 — of, /v0)"] 


, (26) 
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2 respectively. Then hydrodynamic equations (Eq. 16) in the 
3 “Euler step” can be discretized as 
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305 where subscript (i, j) denotes the grid index in the spatial co- 
sos Ordinate system, the superscript n represents the n-th time 
307 Step, T is the time increment, hı and hg are the grid sizes in 
308 the x and y direction, respectively. Note that the mass density 


3% remains constant as shown in Eq. 16, that is pe = pf; 


n n 

ur —u".. 

yi j+1/2 vii—1/2) ut. n 
ER,ij? 


310 C. Lagrange step 


31 In the “Lagrange step”, the fluid is divided into many 
2 pseudo-particles with finite sizes. We assume an internal dis- 
313 tribution function of fluid quantities within a pseudo-particle 
314 [34] 
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aE.) = {(qit1,5-+1 — Gi+1,3)(27 — 69) + qi+1,3} (2E — SE) — { (qi j+1 — qij) (20 — ôn) + qij} (2E — ôE- 1), (29) 


where the pseudo-particle size is normalized, g(£,7) is the 
physical quantity carried by the pseudo-particles°, the coordi- 
nates (£, n) represent the internal position within the pseudo- 
particle , with the origin located at the bottom left-hand cor- 
ner, and (£, ôn) are intervals from the origin of the inter- 
nal coordinate system to the boundaries of the grid in the x 


5 q only represents the velocity u = (ua, Uy ) and specific internal energy e 
of the pseudo-particle, not mass density p, which has been set to constant 
Pp in the initialization. 


32 and y directions, respectively. These are illustrated in Fig. 2, 
323 Where the red lines represent the grids, the red dot represents 
32⁄4 the pseudo-particle center, and the blue square represents the 
32 pseudo-particle size. qp represents the physical quantity car- 
327 ried by the pseudo-particle located at the internal coordinate 
323 center, which can be obtained by integrating q(£, n) from 0 to 
329 1 
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3 
Fig. 2. A distribution of quantities within a pseudo-particle. The red 
dot p represents the pseudo-particle center, the blue square repre- j 


sents the pseudo-particle size, and the black dot is the grid point. 
The label * represents the interval when the pseudo-particle reaches 
a new position. 


qp = ff aenda 
0 


0 
a qi+1,j+1(1 — ôE) (1 — ôn) (30) 
+ qi+1,j(1 — 6€)6n + qi j+1ôE(1 — ôn) + qi jôEðN 
= a(t, 4) 
ZAPI” 


The equations of pseudo-particle motion are 


+ TUp,z, 


(31) 


T TUp,ys 


where (£p, Yp) denote the position of the pseudo-particle cen- 
ter, and (Up x, Up) represent the velocity of the pseudo- 
particle center, as determined by Eq. 30. When the 
pseudo-particle has been advanced and reaches a new po- 
sition as shown in Fig. 2, we recompute the new in- 
tervals from the origin of the internal coordinate, denoted 
as (6€*,d7*). By summing pseudo-particles at new posi- 
tions, new fluid quantities on the grid can be derived. In 
this paper, we present three algorithms, Area-weighting 
method (AWM), Integration-weighting method (IWM) and 
Interpolation-Integration-weighting method (ITWM), for the 
summation process, each of which exhibits different types 
of numerical diffusion [34]. The details of summation algo- 
rithms are shown in Appendix A. The three algorithms have 
different numerical diffusion and noise, and can be selected 
according to the requirements of actual problems. 


D. Parallelization 


The Massage Passing Interface (MPI) is employed in the 
construction of parallelization [35]. The details of parallel 
communication are shown in Appendix B. 
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IV. SIMULATION RESULTS 
A. Far-field X-ray energy deposition 


The module of X-ray energy deposition is crucial for the 
calculation accuracy of the blow-off impulse in Xablation2D 
code. To validate this module, we compared the energy depo- 
sition rate from Xablation2D code with those obtained from 
Monte Carlo code, Geant4 [36-39]. In the Xablation2D sim- 
ulation, the parallel soft X-rays, with an initial energy flux 
po = 418 J/cm?, are incident perpendicularly into a 2D pla- 
nar aluminum (Al) target, with a density p = 2.738 g/cm, 


e and the pulse duration is 50 ns. The energy spectrum of the X- 


rays approximates a black-body spectrum with the radiation 
temperature T = 1 keV. The range of wavelength in energy 


spectrum is AX = 0.1 A to 10 A, discretized into 23 energy 
groups for numerical simulation. The mass absorption coeffi- 
cient u in the Al material for X-rays in each energy group is 
shown in the Tab. 1. 

In the Geant4 code, a parallel black body spectrum photon 
beam with the radiation temperature 1 keV is configured to 
simulate soft X-rays. The interaction between photons and 
materials is simulated by the Livermore low-energy electro- 
magnetic physics model, which is effective within the energy 
range of the photon from 250 eV to 1 GeV. By tracking the 
trajectory of photons within the Al target, it is possible to 
quantify the amount of the deposition energy in the material. 
The energy deposition rates calculated by the two code are 
shown in Fig. 3, where the vertical axis corresponds to the 
energy deposition rate, which signifies the lost ratio of the 
incident energy flux after traversing corresponding distance 
in the material, and the horizontal axis represents the depth 
within the aluminum target, along with the direction of the in- 
cident X-rays. Fig. 3 also shows the relative error of the two 
results by orange dashed line. It is observed that there is the 
largest deviation near the material surface, and then rapidly 
decreases to O(5%) along the incident depth. The reason is 
that the value of the energy deposition rate is small at the ma- 
terial surface, where even a minor deviation on parameters, 
such as the mass absorption coefficients in the database, can 
result in a noticeable relative error. In addition, the depth re- 
quired for 50% X-ray energy to be deposited in the material 
is 4.17 um for Xablation2D and 3.5 wm for Geant4, respec- 
tively, which is consistent with O(5 um) estimated in [19]. 


B. Shear flow simulation 


We use the plane shear flow simulations to verify the nu- 
merical diffusion from the algorithms of AWM, IWM and 
IWM. In the Xablation2D code, the size of the simulation 
domain is Lẹ x Ly = 0.4 x 0.4 cm?, and the domain is 
discretized into 400 x 400 grids. 10 x 10 pseudo-particles 
are allocated for each grid. The whole flow field is divided 
into two layers. The upper fluid is denoted as A and the 
lower fluid as B. The interface between the two fluids is lo- 
cated at y = 0.2 cm. The mass density, velocity in the x 
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TABLE 1. Discretized black-body spectrum from wavelength A = 0.1 A to 10 A and corresponding mass absorption coefficient jz in the Al 


material [19]. 
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Fig. 3. The energy deposition rate varying with the material depth 
along the X-rays incident direction. 


direction and velocity in the y direction of the upper fluid are 
pa = 2.738 ¢/em?, uAx = 3x 10° cm/s and va, = 0 cm/s, 
respectively. The lower fluid has the same mass density as the 
upper pg = pa, a velocity is equal in magnitude but opposite 
in the x direction ugy = —u Ax, and the same velocity in the 
y direction ugy = uay. The simulation employs the periodic 
boundary condition in the x direction, and the open boundary 
condition in the y direction. In addition, the ideal gas EOS is 
adopted, 


p= (y — 1)pe, (32) 


where y = 1.667 is the ratio of specific heat, the initial pres- 
sure p is set to 1 GPa. The fluid velocity distribution in the x 
direction (uz) of shear flow is shown in Fig. 4. 


In Fig. 4, it can be observed that the IWM (or IWM) has 
a better suppression for numerical diffusion than that of the 
AWM. For the AWM in Fig.4 (a), the jumped shear flow ve- 
locity becomes smooth, and this smoothness gradually satu- 
rates into the upper and lower layers due to numerical dif- 
fusion, but for the IWM (or IWM) in Fig.4 (b), the jumped 
velocity can be maintained throughout the entire simulation 
time due to using the integration of the internal distribution 
function in the shear velocity direction. 
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C. Kelvin-Helmholtz instability simulation 


When two contiguous fluids flow with shear velocities, an 
instability can arise at the interface between the fluids as long 
as there is a small fluctuation. This phenomenon is known 
as the Kelvin-Helmholtz instability (KHD [40, 41]. KHI is 
widely observed in natural, including the turbulent mixing of 
fluids in jet streams, the formation of undulates clouds and 
supernova explosions [42—44]. Here, we performed the KHI 
simulations by the Xablation2D code and the open radiation 
magnetohydrodynamics simulation code FLASH4 [45] re- 
spectively, and compared the two simulation results. The ini- 
tial parameters, including the size of the simulation domain, 
the number of grids, and the number of pseudo-particles for 
each grid, are the same as the Section IV B, except for temper- 
ature. The initial material temperature is set to T = 21.6 eV, 
corresponding to the specific internal energy e = 1x 10° J/g. 
The initial pressure determined by EOS is 182.24 GPa. Ini- 
tially, a velocity perturbation in the y direction is introduced 
at the interface of the two fluids as a seed for KHI in the fol- 
lowing form [46] 


u (33) 


y = Uo sin(ka — eee), 


where uo = 1 x 10° cm/s is the initial amplitude, k = 
4r / Ly cm" is the initial wave number. In order to observe 
secondary unstable structures of KHI, it is necessary to re- 
duce numerical diffusion, so we adopt the IWM in the Xab- 
lation2D simulation. In the FLASH4 simulation, we maintain 
consistency in the initialization parameters and perturbation, 
but use the adaptive mesh refinement (AMR). 

Fig. 5 illustrates the mass density the evolution of the lower 
half fluid. The six subfigures in the first row are the results 
obtained from Xablation2D code, and the second row from 
FLASH4 code. The simulation results of the two codes are 
consistent. It can be observed that the disturbance at the fluid 
interface gradually increases, and eventually, the two fluids 
mix and show vortex structures. The IIWM is also employed 
to simulate KHI process with the same parameter in the Xab- 
lation2D. The simulation results are basically consistent with 
FLASH4, except for some secondary structures, as shown in 
Fig. 6. This deviation can be attributed to the relatively higher 
numerical diffusion in the direction of convective velocity in 
the ITWM. 

The mixing width 7 is defined as the difference between 
the maximum position and the initial position of the lower 
half fluid, as shown in Fig. 7 (a). The evolution of the 
mixing width is shown in Fig. 7 (b), where the blue and 
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Fig. 4. The fluid velocity distribution in the x direction in the simulation by two different algorithms at various time intervals. 
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Fig. 5. The mass density distribution of the lower half fluid at various time intervals in the Xablation2D and FLASH4 simulations, respectively 
The first row is corresponding to the result from Xablation2D, and the second is corresponding to the result from FLASH4. 
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Fig. 6. The mass density distribution of the lower half fluid at various time intervals in the Xablation2D simulation by HWM. 
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Fig. 7. The evolution of the mixing width and relative error in the Xablation2D and FLASH4 simulations, respectively. 


FLASH4. The orange dotted lines indicate relative errors be- 
tween two codes, which remain O (5%), indicating a certain 
level of consistency between the two codes. This minor devi- 
ation validates the reliability of the Xablation2D code. It is 
observed that the growth of the mixing width is significantly 
smaller than what is predicted by classical linear theory of the 
KHI. The reason is that the compressibility of the simulation 
and the evolution of the fluid is non-linear due to the strong 
initial perturbation. 


D. Vaporization blow-off impulse 


When X-rays irradiate the solid material, a blow-off im- 
pulse will load on the material surface, as illustrated in Fig. 
8. The energy of the photons in X-rays is primarily absorbed 
by the material through the photoelectric effect and Compton 
scattering. The X-ray energy flux decreases exponentially as 
it penetrates from the surface into the interior of the mate- 
rial. As a result, the specific internal energy of the material at 
the surface is increasing rapidly, which leads to phase transi- 
tion and adiabatic expansion, and finally generates an ablation 
layer. If the deposited energy exceeds the sublimation energy 
of the material, the solid material at the surface changes into 
gas, and forms vaporization zone on the front. The vapor- 
ized material is sprayed into vacuum violently, and generates 
a recoil impulse loading on the remaining material, where the 
recoil impulse is known as the blow-off impulse (BOI). 

In the BOI simulation case, the aluminium (Al) material 
is selected as the target material, the IWM algorithm is em- 
ployed in the Xablation2D code. The expansion and compres- 
sion of material are described by PUFF EOS (Eq. 10) and 
Giineisen EOS (Eq. 12), as discussed in Section IC. The 
physical property parameters of the Al material in the sim- 
ulation are shown in the Tab. 2. The simulation domain has a 
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Fig. 8. A simple illustration of the far-field X-ray ablation. The 
shaded area represents the blow-off impulse, and the black solid 
curve is the energy deposition curve. 


sos Size of Ly x Ly = 0.80.8 cm?, and is divided into 400 x 400 
so7 grids. We allocate 100 x 100 pseudo-particles for each grid 
sos near the surface of the target material, while 5 x 5 pseudo- 
so particles in other region.Two different geometric configura- 
510 tions of the target are simulated, one is a semi-infinite slab 
s1 and the other is a cylinder. 
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po (g/cm?) [2.738 || Yo (GPa) | 0.7 |les (kJ/g)| 10.89 
Co (mm/s) |5.328]| G (GPa) | 27 y 1.667 
À 1.338 To l2as|| N  J1.265 


TABLE 2. The physical property parameters of the Al material 
in the simulation. po is the mass density at at room temperature 
and atmospheric pressure, G is the shear modulus, Yo is the yield 
strength, es is the sublimation energy, y is the specific heat ratio, T'o 
is the Griineisen coefficient, A and Co are Hugoniot parameters, and 
N =CG/Toes . 


1. Semi-Infinite Slab 


The semi-Infinite slab is located in the region where x > 
0.4 cm, while the region where x < 0.4 cm is filled with low- 
density gas to maintain the numerical stability of the code. 
The periodic boundary condition is applied in the y direction, 
and the open boundary condition is applied in the x direction. 
The parallel X-rays have a radiation temperature T = 1 keV, 
and incident on the slab along the normal direction. The en- 


BBAY : 


Whitener : 


MBBAY : 


where e(x) is energy deposition profile, em represents the 
melting energy, a is a correction parameter ranging from 1 
to v2, and zo is the thickness of the sublimation layer de- 
termined by setting e(x) = es. In the simulation, we set 
em = 3.975 kJ/g and a = 1.1, respectively. It should be 
noted that the BOI models only estimate the instantaneous 
deposition energy for irradiation, resulting in horizontal dot- 
ted lines in Fig. 10. While the BOIs calculated by the Xabla- 
tion2D and 1D Lagrangian codes, shown as the blue dashed 
and magenta solid lines in Fig. 10, has a growth time of 
O(0.1 us), it demonstrates a tendency to stabilize, which cor- 
responds to the characteristic time in which the material com- 
pletely sublimates at the surface. It is also observed that the 
results obtained from the BBAY model and MBBAY model 
align well with the stable values of the numerical simulations. 


Furthermore, we compared the simulation results of the 
Xablation2D code with the published experimental results to 
verify the reliability of the code. The Ref. [48] provides three 
measurements of the BOI produced by X-rays irradiating a 
flat Al material. The X-ray parameters and measured impulse 
values in experiments are presented in the Tab. 3. We simu- 
late these three experiments by Xablation2D with the experi- 
mental parameters. The BOIs of the simulations are shown in 
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ergy spectrum is discretized into 23 energy groups as in Tab. 
1. The pulse duration and initial energy flux are 50 ns and 
418 J/ cm”, respectively. Fig. 9 illustrates the evolution of the 
mass density and pressure. It is observed that the material sur- 
face is vaporized by the X-rays irradiation, and then expand 
outward into surroundings. The material expanding generates 
a blow-off impulse (BOI), and gives rise to a thermal shock 
wave propagating inward. The accumulation of density in the 
left region is a consequence of the presence of low-density 
gas in the initialization. 

For the semi-infinite slab configuration, the impulse distri- 
bution is uniform in the y direction, it is possible to obtain 
the evolution of one-dimensional average BOI by integrating 
it along the y direction, and dividing it by the length Ly. This 
result shown as the blue dashed line in Fig. 10. To verify its 
correctness, we develop a one-dimensional Lagrangian code 
referred to Ref. [47] to calculate the BOI under the same pa- 
rameters. The result is shown as the magenta solid line, which 
is consistent with the Xablation2D code. In addition, we also 
compare the simulation results with three BOI models shown 
as the dotted lines in Fig. 10. Their analytical formulas are 
written as follows [9], 


1/2 
(34) 


1/2 
i) pax a} ; 


Em 


Fig. 11 (a). The stable values of these BOIs are 97.84, 121.82, 
125.36 Pa-s, respectively. Fig. 11 (b) shows the comparison 
between the simulation results and the measured BOIs. It is 
observed that there are two consistent BOIs, and the relative 
error between the simulation and the experimental results are 
negligible in case of No.01154 and No.01170. Although the 
relative error in case of No.01171 is larger, it is still within 


O(20%). 


Number 01154/01170/01171 

Radiation temperature T' (keV) |] 0.21 | 0.227 | 0.211 
Pulse duration 7o (ns) 53 36 44 
Initial flux ®o (J/cm?) 163 | 181 | 192 

Measuring impulse 7 (Pa-s) || 99.1 | 118.5 | 162.2 


TABLE 3. The X-ray parameters and measuring impulse in experi- 
ments. 
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Fig. 9. The mass density and pressure distribution in the semi-inifnite slab at various time intervals. 
ss7 thermal shock wave moving toward the center of the cylinder. 
ssa When the parallel X-ray irradiates on a curved surface, the 
sæ density of the deposited energy at the material surface is non- 
QO yor 8 E ee arene erie er SvaPTeSer so uniform. The BOI should be a function of the radial direction 
s92 Of the cylinder. Refs. [10, 11, 21, 49] propose that if the en- 
sə ergy deposition density has a cosine profile on the cylinder 
3004 sə surface, the variation of the BOI is proportional to the cosine 
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Fig. 10. The blow-off impulse varying with the simulation time. 


2. Cylinder 


In the cylinder configuration, the target is positioned in the 
center of the simulation domain. The central coordinates of 
the cylinder is at x = y = 0.4 cm, and the radius of the target 
is 0.2 cm. The remaining areas of the simulation domain are 
filled with a low-density gas. The physical property parame- 
ters and X-ray parameters used in this configuration are iden- 
tical to those employed in the semi-infinite slab configuration. 
The mass density and pressure distribution in the simulation 
are shown in Fig. 12. The X-ray irradiates the left side of the 
cylinder. The vaporized material is produced on the left sur- 
face, ejects violently, and generates the BOI, which drives a 
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of the polar angle 


I = Ío: cos, (35) 
where @ is the polar angle defined in Fig. 13 (a), Io represents 
the BOI at the 0 = 0°. 

Fig. 13 (b), (c), and (d) show the distribution of the BOI 
via the polar angle at different time intervals. The BOI pro- 
files basically satisfies the cosine law, but small deviations 
are also exist. The main reason is that the thickness of X- 
ray energy deposition cannot be completely ignored. The 
deviation between the BOI curve and the cosine function in- 
creases as the polar angle increases from 0° to 90°, especially 
at 9 = 90°, where specific internal energy of the material is 
not zero due to the penetration of X-rays into the material, so 
the BOI is also not zero. The simulation results confirm this 
phenomenon. 


V. CONCLUSION 


We develop a novel parallel two-dimensional PIC code, 
Xablation2D, for calculating the BOI of far-field X-ray ir- 
radiating materials. The code is mainly composed of three 
modules: energy deposition module, EOS module and ideal 
hydrodynamics module. In ideal hydrodynamics module, 
the solution of hydrodynamics equations is divided into two 
steps: “Euler step” and “Lagrange step” according to the 


sts Harlow’s splitting scheme. The introduced pseudo-particle 
si9 is responsible for compensating for the transport effect. We 631 
620 present three new summation algorithms, AWM, IWM and 
e21 IIWM, to map the physical quantity carried by the pseudo- 
e22 particle into the grid. In contrast to the conventional finite 
difference or finite element methods, the new PIC method sig- 
nificantly reduces the dependence on grid shape and is well- 
suited for the calculation of large deformation problems in 637 
materials with complex configurations. In order to verify the 
reliability of the Xablation2D code, we use some open-source 
codes for comparison. The soft X-ray energy deposition rate 
is simulated by Geant4, and the shear flow and the KHI prob- 
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(b) The simulation results vs the experiment results 


Fig. 11. Comparison between the Xablation2D simulations and experiments. 
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Fig. 12. The mass density and pressure distribution in the cylinder at various time intervals. 
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by FLASH4, respectively. The simulation 


results exhibit a minor discrepancy compared to our code. Fi- 


the BOI in two geometric configurations 
under far-field X-ray irradiation. The re- 


sults are well consistent with experimental, analytical, and 


other simulation results. 
shock wave propagates inside the material caused by X-ray 


It is also observed that a thermal 


increased energy flux, the overall process 


becomes more intricate. The material will be highly ionized, 
e4 and the energy deposition process involves the interaction be- 
tween the radiation field and the plasma. In the future work, 
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we plan to develop a radiation transport module in the Xabla- ess uct of qp and area as follows, 


tion2D code for the simulation of ultra-intense X-ray irradi- 


ating materials. 
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Appendix A: Summation algorithms 


(pq)i+1,j+1 = 


pE 


pE 


i,j) 

XO ppd (1 — 6€")6n", 
i+1,7) 

NO ppdpdé* (1 ôn"), 
4,j+1) 


So PpGp(1 — 6&*)(1 ôn"). 


i+1,j+1) 


(Al) 


660 Where p € (i, j) denotes the pseudo-particle whose cen- 
661 ter is located in the (i, 7) grid. The AWM is a straightfor- 
662 ward approach and has high computational efficiency, but 
663 introduce zero order numerical diffusion effects. 
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4 (2) Integration-weighting method (IWM): The IWM main- 


665 tains the internal distribution function q(£, n) during the 


666 motion of the pseudo-particle. The new physical quanti- 

(1) Area-weighting method (AWM): This approach as- 67 ties on the grid can be obtained by integrating the func- 
sumes that the distribution function of the physical quan- sss tion q(€,7) over each new area. This area is defined by 
tity g(&,7) carried by the pseudo-particle collapses into «6s the overlap of the pseudo-particle size and the grid. For 
the value located at the internal coordinate center, that is 670 instance, the area of overlap between the pseudo-particle 

qp. Therefore, the new fluid physical quantities on the 67 size and the (i, j) grid can be represented by 6&* x ôn* 


grid can be simply obtained by a summation of the prod- 672 as shown in Fig. 2. As a result, the new fluid physical 
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pe (i,j) 
1 ôn 
(is = X m» f q(, n)d&dn, 
i415 Teg 
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pe (i,j +1) 0 én* 
1 1 
(einige = >» F, | q(é,n)dédn. 
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The IWM has the second-order numerical diffusion, 
which is much better than that of the AWM. However, 
this approach would introduce some non-physical numer- 
ical noise, causing the fluid motion become somewhat un- 
stable. 


(3) Interpolation-Integration-weighting method (ITWM): 
The ITWM is a combination of the AWM (used in the con- 
vective velocity direction) and IWM (used in the shear 
velocity direction). For instance, if the g(, n) represents 
the fluid velocity in the x direction (uz), then an inverse 
linear interpolation is employed to q(€,7) in the x direc- 
tion, and followed by the integration of q(€,7) in the y 
direction. The specific procedure are as follows. Initially, 
we integrates the q(€, 7) over € from 0 to 1, 


a(n) = J aE, nde 


0 
= [dea gaa (1 — SE) + qi j+1ôE] 
— [qi+1,;(1 — 8E) + qi jôE]} (2n — ôn) 
+ qi+1,;(1 — 6€) + qi jôÈ 


= (57). 


(A3) 


Then, we integrate q(7)) over the length where the overlap 
of the pseudo-particle size and the grid in the y direction, 
as shown in Fig. 2, to obtain two new intermediate phys- 
ical quantities 
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Finally, new physical quantities are determined by 


(Paji; = XO pride", 


pe (i,j) 
(eq)it1,5 = 5 PpG (1 — 6&*), 

pe(i+1,j) (A5) 
(dij = XO Pit”, 

pE(ij+1) 


(pq)i+1,j+1 = Pp (1 — He"). 


D 


pE(i+1,j+1) 


Similarly, we can also implement this approach for the 
fluid velocity in the y direction (uy). 


Appendix B: Parallel communication 


MPI adopts an explicit message passing architecture, ne- 
cessitating the explicit sending and receiving of messages 
to facilitate data exchange among processors.Each individual 
parallel process possesses its own memory space. Processors 
access the memory space of each other through explicit mes- 
sage passing. This design exhibits robust parallelism and is 
particularly well-suited for the implementation of large-scale 
parallel algorithms. 


Fig. 14. The parallel division of the simulation region. A total of 
N = 12 processes, numbered from 0 to 11, are assigned to the 
simulation regions. The number of processes assigned to the grid 
along the x and y direction are 4 and 3, respectively. 


In our code, the simulation region is divided into N parts 
based on a grid structure, and each part is assigned to a indi- 
vidual process for calculation. The total number of processes 
N = ppx x ppy, where ppx and ppy are the number of pro- 
cesses assigned to the grid along the x and y direction, re- 
spectively, as shown in Fig. 14. A layer of grid, named as 
guard grid or ghost grid, is set up at the junction of adjacent 
parallel regions to store the physical quantities received from 
the adjacent processes. The exchange of physical quantities 
on the guard grids is achieved through two MPI communica- 
tion subroutine: MPI_SEND and MPI_RECV, as illustrated 
in Fig. 15. 
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Fig. 15. MPI is utilized to update the guard grids. The dashed boxes 
represent parallel regions assigned to processors. The black and grey 
solid circles represent the grid points within the parallel region and 
guard points surrounding the parallel region, respectively. The red 
and blue arrow lines depict the communication that occurs between 
internal grid points and guard points across adjacent processors. 


In the initialization, all pseudo-particles are placed within 
the parallel region and guard grids. When pseudo-particles 
are pushed out of the parallel region, they will be trans- 
ferred to another parallel region by MPI communication. This 
communication can be realized between any two processors. 
Firstly, all communicating pseudo-particles are gathered into 
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